Staggered-grid Finite-difference Scheme Research Articles

A stable staggered-grid finite-difference scheme for acoustic modeling beyond conventional stability limit

Staggered-grid finite-difference (SGFD) schemes have been widely used in acoustic wave modeling for geophysical problems. Many improved methods are proposed to enhance the accuracy of numerical modeling. However, these methods are inevitably limited by the maximum Courant-Friedrichs-Lewy (CFL) numbers, making them unstable when modeling with large time sampling intervals or small grid spacings. To solve this problem, we extend a stable SGFD scheme by controlling SGFD dispersion relations and maximizing the maximum CFL numbers. First, to improve modeling stability, we minimize the error between the FD dispersion relation and the exact relation in the given wave-number region, and make the FD dispersion approach a given function outside the given wave-number area, thus breaking the conventional limits of the maximum CFL number. Second, to obtain high modeling accuracy, we use the SGFD scheme based on the Remez algorithm to compute the FD coefficients. In addition, the hybrid absorbing boundary condition is adopted to suppress boundary reflections and we find a suitable weighting coefficient for the proposed scheme. Theoretical derivation and numerical modeling demonstrate that the proposed scheme can maintain high accuracy in the modeling process and the value of the maximum CFL number of the proposed scheme can exceed that of the conventional SGFD scheme when adopting a small maximum effective wavenumber, indicating that the proposed scheme improves stability during the modeling.

Mimetic finite differences for boundaries misaligned with grid nodes

This paper considers the problem of mimetic staggered-grid finite difference scheme construction in the case of a boundary misaligned with the grid points. We show that depending on the reciprocal positions of the nearest grid point and whether the point is inside or outside the domain, two cases should be considered separately. In both cases, we present general rules to construct the mimetic difference operators to approximate the gradient and divergence operators. In addition, we investigate the stability criterion of the explicit in time scheme, which employs derived operators to approximate spatial derivatives. We show that, depending on the reciprocal positions of the nearest grid point, the Courant number may be smaller by 7% than that of the initial value problem approximation.

Free-surface elastic wave injection and wavefield reconstruction with applications to elastic reverse time migration

Accurately injecting and extrapolating multicomponent seismic data is a major step in elastic reverse time migration (RTM) and full-waveform inversion. The theoretical basis of this process is the seismic representation theorem. We propose a numerical approach based on the collocated-grid MacCormack finite-difference (FD) scheme to implement the injection and extrapolation. The MacCormack FD scheme defines wavefield parameters on the same grid. When dealing with the boundary condition in the representation theorem, it is more advantageous than the commonly used staggered-grid FD scheme which defines the stress and velocity on offset grids; thus, interpolations are required to link the boundary values to nearby nodes. In contrast, with the MacCormack scheme, the traction and particle velocity data can be directly injected into the boundary nodes as equivalent concentrated forces, force dipoles, and double couples. The accuracy of this approach is validated by comparing the regenerated wavefield with the original wavefield. Our primary purpose is to develop an injection method for elastic RTM, but this approach also can be used to numerically investigate other interesting issues in the elastic wave injection, including errors and artifacts caused by injection from unclosed boundaries, incomplete boundary values, and uncertain velocity models. In particular, seismic data are usually acquired on the ground surface where the traction-free condition is naturally maintained. We simulate a land acquisition scenario by generating a synthetic multicomponent data set on a traction-free surface. The reconstructed wavefields from velocity data are compared with the original forward-propagated waves to evaluate the accuracy of the boundary treatments. Two multicomponent elastic RTM examples are presented to demonstrate the applications of the proposed injection method. By combining the regenerated receiver-side wavefield and the vectorized imaging condition, the PP, PS, SP, and SS images are calculated.

Controlled-source electromagnetic modeling using a high-order finite-difference time-domain method on a nonuniform grid

Simulation of 3D low-frequency electromagnetic (EM) fields propagating in the earth is computationally expensive. We have developed a fictitious wave-domain high-order finite-difference time-domain (FDTD) modeling method on nonuniform grids to compute frequency-domain 3D controlled-source EM data. The method overcomes the inconsistency issue widely present in the conventional second-order staggered-grid finite-difference scheme over a nonuniform grid, achieving high accuracy with an arbitrarily high-order scheme. The finite-difference coefficients adaptive to the node spacings can be accurately computed by inverting a Vandermonde matrix system using an efficient algorithm. A generic stability condition applicable to nonuniform grids is established, revealing the dependence of the time step and these finite-difference coefficients. A recursion scheme using fixed-point iterations is designed to determine the stretching factor to generate the optimal nonuniform grid. The grid stretching in our method reduces the number of grid points required in the discretization, making it more efficient than the standard high-order FDTD with a densely sampled uniform grid. Instead of stretching in vertical and horizontal directions, better accuracy of our method is observed when the grid is stretched along the depth without horizontal stretching. The efficiency and accuracy of our method are demonstrated by numerical examples.

Three-dimensional elastic reverse-time migration using a high-order temporal and spatial staggered-grid finite-difference scheme

Three-dimensional (3D) elastic reverse-time migration (ERTM) can image the subsurface 3D seismic structures, and it is an important tool for the Earth’s interior imaging. A common simulation kernel used in 3D ERTM is the current staggered-grid finite-difference (SGFD) method of the first-order elastic wave equation. However, the mere second-order accuracy in time of the current SGFD method can bring non-negligible time dispersion, which reduces the simulation accuracy and further leads to the distortion of the imaging results. This paper proposes a vector-based 3D ERTM using the high-order accuracy SGFD method in time to obtain high-accuracy images. This approach is a new high-resolution ERTM workflow that improves the imaging accuracy of conventional ERTM from numerical simulation. The proposed ERTM workflow is established on a quasi-stress–velocity wave equation and its vector wavefield decomposition form. Advanced SGFD schemes and their corresponding coefficients with fourth-order temporal accuracy solve the quasi-linear wave equation system. The normalized dot product imaging condition produces high-quality images using high-accuracy vector wavefields solved using the SGFD method. Through the numerical examples, we test the simulation efficiency and analyze how temporal accuracy in numerical simulations affects migration imaging quality. We include that the proposed method obtains highly accurate images.

Non-artifact vector P- and S-wave separation for elastic reverse time migration

Elastic reverse time migration (RTM) uses the elastic wave equation to extrapolate multicomponent seismic data to the subsurface and separate the elastic wavefield into P- and S-waves. P- and S-wave separation is a necessary step in elastic RTM to avoid crosstalk between coupled wavefields. However, the current curl-divergence operator-based separation method has a polarity reversal problem in PS imaging, and vector separation methods often have separation artifacts at the interface, which affects the quality of the imaging stack. We propose a non-artifact P- and S-wave separation method based on the first-order velocity-strain equation. This equation is used for wavefield extrapolation and separation in the first-order staggered-grid finite-difference scheme, and the storage and calculation amounts are consistent with the classical first-order velocity-stress equation. The separation equation does not calculate the partial derivatives of the elastic parameters, and thus, there is no artifact in the separated P- and S-waves. During wavefield extrapolation, the dynamic characteristics of the reflected wave undergo some changes, but the transmitted wavefield is accurate; therefore, it does not affect the dynamic characteristics of the final migration imaging. Through numerical examples of 2D simple models, part SEAM model, BP model, and 3D 4-layer model, different wavefield separation methods and corresponding elastic RTM imaging results are analyzed. We found that the velocity-strain based elastic RTM can image subsurface structures well, without spike artifacts caused by separation artifacts, and without polarity reversal phenomenon of the PS imaging.

Wave Equation Numerical Simulation and RTM With Mixed Staggered-Grid Finite-Difference Schemes

For the conventional staggered-grid finite-difference scheme (C-SFD), although the spatial finite-difference (FD) operator can reach 2Mth-order accuracy, the FD discrete wave equation is the only second-order accuracy, leading to low modeling accuracy and poor stability. We proposed a new mixed staggered-grid finite-difference scheme (M-SFD) by constructing the spatial FD operator using axial and off-axial grid points jointly to approximate the first-order spatial partial derivative. This scheme is suitable for modeling the stress–velocity acoustic and elastic wave equation. Then, based on the time–space domain dispersion relation and the Taylor series expansion, we derived the analytical expression of the FD coefficients. Theoretically, the FD discrete acoustic wave equation and P- or S-wave in the FD discrete elastic wave equation given by M-SFD can reach the arbitrary even-order accuracy. For acoustic wave modeling, with almost identical computational costs, M-SFD can achieve higher modeling accuracy than C-SFD. Moreover, with a larger time step used in M-SFD than that used in C-SFD, M-SFD can achieve higher computational efficiency and reach higher modeling accuracy. For elastic wave simulation, compared to C-SFD, M-SFD can obtain higher modeling accuracy with almost the same computational efficiency when the FD coefficients are calculated based on the S-wave time–space domain dispersion relation. Solving the split elastic wave equation with M-SFD can further improve the modeling accuracy but will decrease the efficiency and increase the memory usage as well. Stability analysis shows that M-SFD has better stability than C-SFD for both acoustic and elastic wave simulations. Applying M-SFD to reverse time migration (RTM), the imaging artifacts caused by the numerical dispersion are effectively eliminated, which improves the imaging accuracy and resolution of deep formation.

A hybrid explicit implicit staggered grid finite-difference scheme for the first-order acoustic wave equation modeling

Implicit staggered-grid finite-difference (SGFD) methods are widely used for the first-order acoustic wave-equation modeling. The identical implicit SGFD operator is commonly used for all of the first-order spatial derivatives in the first-order acoustic wave-equation. In this paper, we propose a hybrid explicit implicit SGFD (HEI-SGFD) scheme which could simultaneously preserve the wave-equation simulation accuracy and increase the wave-equation simulation speed. We use a second-order explicit SGFD operator for half of the first-order spatial derivatives in the first-order acoustic wave-equation. At the same time, we use the implicit SGFD operator with added points in the diagonal direction for the other first-order spatial derivatives in the first-order acoustic wave-equation. The proposed HEI-SGFD scheme nearly doubles the wave-equation simulation speed compared to the implicit SGFD schemes. In essence, the proposed HEI-SGFD scheme is equivalent to the second-order FD scheme with ordinary grid format. We then determine the HEI-SGFD coefficients in the time–space domain by minimizing the phase velocity error using the least-squares method. Finally, the effectiveness of the proposed method is demonstrated by dispersion analysis and numerical simulations.

A novel equivalent staggered-grid finite-difference scheme and its optimization strategy for variable-density acoustic wave modelling

Compared with the standard staggered-grid finite-difference (FD) methods, equivalent staggered-grid (ESG) ones can significantly reduce the computational memory for acoustic wave modelling in the variable-density media. To further enhance the simulation efficiency and accuracy, one way is to optimize the FD coefficients, another way is to design new FD stencils. In this paper, we propose a modified ESG (M-ESG) scheme which can significantly accelerate the wavefield simulation process while preserving or even improving the modelling accuracy. We calculate the FD coefficients by approximating the temporal and spatial derivatives simultaneously based on time–space domain (TS-D) dispersion relation of the discrete wave equation. Our M-ESG scheme in the TS-D can maintain basically the same accuracy as the conventional ESG (C-ESG) one when the FD coefficients are derived by the Taylor-series expansion (TE) approach. Note that the TS-D dispersion relation is nonlinear with respect to the FD coefficients of the C-ESG scheme, so it is difficult to obtain the optimized FD coefficients for the discrete wave equation. However, we can minimize the L2-norm error of the dispersion relation based on our M-ESG scheme to implement a linear FD coefficients optimization strategy, which is easy and efficient. Comparisons with TE- and optimization-based C-ESG schemes demonstrate the accuracy, stability, and efficiency superiorities of our TE- and optimization-based M-ESG ones.

Time–space domain scalar wave modeling by a novel hybrid staggered-grid finite-difference method with high temporal and spatial accuracies

Staggered-grid finite-difference (SFD) scheme is favored in wave equation simulation due to its superior accuracy and stability to center-grid finite-difference (CFD) scheme. However, for the scalar wave equation (SWE) modeling, the conventional SFD (CSFD) scheme only reaches second-order accuracy in space and time for the discrete SWE; the recently developed modified SFD (MSFD) scheme improves the accuracy to (2N)th-order (N<4) but is costly because the MSFD scheme is coined by adding many extra grid points to the CSFD scheme. To tackle these issues, we develop a cost-effective hybrid SFD (HSFD) scheme, which combines the features of the CSFD and MSFD schemes; we prove that the new HSFD scheme can simultaneously reach (2N)th-order accuracy in space and time for the discrete wave equation. In addition, to deal with the optimization difficulties due to the nonlinear dispersion relation of the SFD schemes, we propose a two-step linear optimization method to improve the accuracy of the new HSFD scheme. The analyses on dispersion, stability properties and numerical simulation examples demonstrate that the new HSFD scheme owns better accuracy and stability than the CSFD scheme. Computational cost analysis shows that the HSFD scheme can be more efficient than the MSFD scheme because it only requires half the additional grid points.

Propagation behavior of Love-type waves in poro-elastic medium using staggered grid finite difference Scheme

The analysis of the propagation behaviour of seismic waves in porous medium requires a mathematical backup. That is pictured here. The present study investigates the propagation behavior of Love-type waves in poro-elastic medium. A staggered grid finite-difference (SFD) scheme in time-space domain formulation for the Biot’s equation of poro-elasticity, is presented and applied to the problem. Complete theoretical seismograms for the horizontal and vertical components of displacement are obtained. The dispersion curves are evaluated considering the material parameters of different models. The stability analysis of the above numerical scheme is deliberated. From different graphs, it is noticed that there is difference in phase velocity for various models. By observing the curves behaviour, we can understand the nature of the composite structure. Our outcomes also endorse the finite difference modeling is an important numerical tool to acquire knowledge about the transmission of seismic waves in porous medium. The agreement is excellent.

Acoustic wave propagation with new spatial implicit and temporal high-order staggered-grid finite-difference schemes

AbstractThe implicit staggered-grid (SG) finite-difference (FD) method can obtain significant improvement in spatial accuracy for performing numerical simulations of wave equations. Normally, the second-order central grid FD formulas are used to approximate the temporal derivatives, and a relatively fine time step has to be used to reduce the temporal dispersion. To obtain high accuracy both in space and time, we propose a new spatial implicit and temporal high-order SG FD stencil in the time–space domain by incorporating some additional grid points to the conventional implicit FD one. Instead of attaining the implicit FD coefficients by approximating spatial derivatives only, we calculate the coefficients by approximating the temporal and spatial derivatives simultaneously through matching the dispersion formula of the seismic wave equation and compute the FD coefficients of our new stencil by two schemes. The first one is adopting a variable substitution-based Taylor-series expansion (TE) to derive the FD coefficients, which can attain (2M + 2)th-order spatial accuracy and (2N)th-order temporal accuracy. Note that the dispersion formula of our new stencil is non-linear with respect to the axial and off-axial FD coefficients, it is complicated to obtain the optimal spatial and temporal FD coefficients simultaneously. To tackle the issue, we further develop a linear optimisation strategy by minimising the L2-norm errors of the dispersion formula to further improve the accuracy. Dispersion analysis, stability analysis and modelling examples demonstrate the accuracy, stability and efficiency advantages of our two new schemes.

Reverse Time Migration Based on the Pseudo-Space-Domain First-Order Velocity-Stress Acoustic Wave Equation

Reverse time migration (RTM) is an ideal seismic imaging method for complex structures. However, in conventional RTM based on rectangular mesh discretization, the medium interfaces are usually distorted. Besides, reflected waves generated by the two-way wave equation can cause artifacts during imaging. To overcome these problems, a high-order finite-difference (FD) scheme and stability condition for the pseudo-space-domain first-order velocity-stress acoustic wave equation were derived, and based on the staggered-grid FD scheme, the RTM of the pseudo-space-domain acoustic wave equation was implemented. Model experiments showed that the proposed RTM of the pseudo-space-domain acoustic wave equation could systematically avoid the interface distortion problem when the velocity interfaces were considered to compute the pseudo-space-domain intervals. Moreover, this method could effectively suppress the false scattering of dipping interfaces and reflections during wavefield extrapolation, thereby reducing migration artifacts on the profile and significantly improving the quality of migration imaging.

First-order particle velocity equations of decoupled P- and S-wavefields and their application in elastic reverse time migration

The separation of P- and S-wavefields is considered to be an effective approach for eliminating wave-mode crosstalk in elastic reverse time migration (RTM). At present, the Helmholtz decomposition method is widely used for isotropic media. However, it tends to change the amplitudes and phases of the separated wavefields compared with the original wavefields. Other methods used to obtain pure P- and S-wavefields include the application of the elastic wave equations of the decoupled wavefields. To achieve high computational accuracy, staggered-grid finite-difference (FD) schemes are usually used to numerically solve the equations by introducing an additional stress variable. However, the computational cost of this method is high because a conventional hybrid wavefield (the P- and S-wavefields are mixed together) simulation must be created before the P- and S-wavefields can be calculated. We have developed the first-order particle velocity equations to reduce the computational cost. The equations can describe four types of particle-velocity wavefields: the vector P-wavefield, the scalar P-wavefield, the vector S-wavefield, and the vector S-wavefield rotated in the direction of the curl factor. Without introducing the stress variable, only the four types of particle velocity variables are used to construct the staggered-grid FD schemes; therefore, the computational cost is reduced. We also develop an algorithm to calculate the P and S propagation vectors using the four particle velocities, which is simpler than the Poynting vector. Finally, we apply the velocity equations and propagation vectors to elastic RTM and angle-domain common-image gather computations. These numerical examples illustrate the efficiency of our methods.

A comparative study on the stress image and adaptive parameter-modified methods for implementing free surface boundary conditions in elastic wave numerical modeling

The stress image method is a simple and effective approach for implementing the free surface boundary conditions in elastic wave numerical modeling. This method assumes that the normal and shear stresses perpendicular to the free surface are antisymmetric with respect to the free surface. In this way, the values of the normal and shear stresses above the free surface can be updated. However, the stress image method is based on an intuitive viewpoint and lacks a physical foundation. The adaptive parameter-modified method is a recently proposed approach for implementing the free surface boundary conditions. Through adaptive modifications of the density and elastic parameters at the free surface, the implementation of the free surface boundary conditions is achieved. Based on the adaptive parameter-modified method, a new interpretation of the stress image method is used. We determine that the stress image method is equivalent to the adaptive parameter-modified method in terms of the staggered-grid finite-difference scheme for the elastic wave equation in displacement form. This result provides a physical foundation and explanation of the stress image method. Therefore, we can further develop the stress image method from a physical viewpoint. Numerical examples are also developed to perform the theoretical analysis.

Analysis of Seismic Wavefield Characteristics in 3D Tunnel Models Based on the 3D Staggered-Grid Finite-Difference Scheme in the Cylindrical Coordinate System

The accurate prediction of the geological conditions ahead of a tunnel plays an important role in tunnel construction. Among all forward geological prospecting methods, the seismic detection method is widely applied. However, due to the characteristics of the tunnel and the complexity of the geological conditions, the seismic wavefield is complicated. Carrying out a more realistic forward modeling method is vital for fully understanding the law of seismic wave propagation and the characteristics of seismic wavefield in the tunnel. In this paper, the 3D staggered-grid finite-difference scheme in the cylindrical coordinate system based on the decoupled nonconversion elastic wave equation is used to carry out the numerical simulation. This method can avoid the diffraction interferences produced at the edges of the tunnel face in the Cartesian coordinate system. Based on this forward modeling method, the characteristics of wavefield and propagation laws of seismic waves under three kinds of common typical unfavorable geological models were explored, which can provide theoretical guidance to seismic data interpretation of tunnel seismic forward prospecting in practice.

Elastic Wave Modeling With High-Order Temporal and Spatial Accuracies by a Selectively Modified and Linearly Optimized Staggered-Grid Finite-Difference Scheme

High-order staggered-grid finite-difference (SFD) schemes are preferred in elastic wave simulation for geophysical problems because they decrease the accumulation of error from grid dispersion. However, most SFD approaches reach high-order spatial but limited temporal accuracy. To tackle the issue, we develop a novel temporal and spatial high-accuracy elastic SFD scheme by selectively modifying the spatial operators of the original SFD stencil. This modification has three main advantages. First, it facilitates the design of a new SFD stencil with temporal and spatial accuracies to arbitrary even-order by a Taylor-series expansion method. Second, it helps boost the accuracy further by implementing a linear optimization method. Third, the new selectively modified SFD (SMSFD) stencil needs fewer float-point operations (FPOs) than the existing temporal high-order SFD stencil. We compare our new SMSFD scheme with spatial high-order and temporal–spatial high-order SFD schemes and show that our new elastic SMSFD scheme possesses better accuracy and stability and requires fewer FPOs than these methods.

3D general-measure inversion of crosswell EM data using a direct solver

Abstract We present a three-dimensional (3D) general-measure inversion scheme of crosswell electromagnetic (EM) data in the frequency domain with a direct forward solver. In the forward problem, we discretised the EM Helmholtz equation by the staggered-grid finite difference (SGFD) scheme and solved it using the Intel MKL PARDISO direct solver. By applying a direct solver, we simultaneously solved the multisource forward problems at a given frequency. In the inversion, we integrated a general measure of data misfit and model constraints with linearised least-squares inversion. We reconstructed a model with blocky features by selecting the appropriate measure parameters and model constraints. We used the adjoint equation method to explicitly calculate the Jacobian matrix, which facilitated the determination of an appropriate initial value for the regularisation coefficient in the objective function. We illustrated the inversion scheme using synthetic crosswell EM data with a general measure, the L2 norm, and, specifically, two mixed norms.

3-D Crosswell Electromagnetic Inversion Based on General Measures

The crosshole electromagnetic (EM) method is used in reservoir development monitoring and evaluation. It reflects the resistivity distribution between wells in the form of EM response changes and describes the spatial distribution of pore fluids in subsurface reservoirs based on the correlation between resistivity and pore fluids. This technique has been used to describe reservoirs, monitor steam and brine flooding to improve oil recovery, optimize the design of infilled wells, and verify the effectiveness of reservoir numerical simulations. However, due to the complexity of the electrical distribution of underground reservoirs and the limitations of the crosswell EM response range, there is still a lack of practical and effective crosswell EM imaging methods. This article presents a three-dimensional EM inversion scheme for interpreting crosswell data in the frequency domain with general measures. For the forward modeling, the equation for the secondary electric field was discretized using the staggered-grid finite difference scheme and was solved using the Krylov subspace iterative method. For the inversion problem, the general measure of the data misfit and regularization constraints were integrated using the Gauss–Newton minimization method to form an objective function. By selecting appropriate parameters and regularization constraints, the sharp boundaries and blocky features of the inverted target were recovered. The sensitivity was calculated using the reciprocity relationship, and the method was found to be beneficial to determining appropriate initial values of regularization coefficients. Two synthetic data sets and a field data inversion of a crosswell configuration demonstrate the feasibility of this scheme.

An improved rotated staggered grid finite difference scheme in coal seam

Coal seam is thin compared to the wavelength of seismic wave and usually shows strong anisotropy. It may form special geological bodies such as goaf and collapse in coal mining. The existence of these phenomena may lead to instability in numerical simulation of goaf area in coal seam. And the calculation speed of simulation is always a factor that restricting the development of simulation technique. To improve the accuracy and efficiency of seismic numerical simulation of goaf area, an improved vacuum method has been imported in a rotated staggered grid scheme and the calculation is implemented by combining parallel computing and task parallelism. It makes sure that the numerical simulation method can be utilized in a geological model with large differences in elastic parameters among layers and improve the performance of a parallel application by enabling full use of processor resources to expedite calculation speed. We set up anisotropic coal seam models, and then analyze the characteristics of synthetic seismograms and snapshots of different goaf areas with or without collapse by modeling. Results show that the proposed method can accurately simulate the goaf area and the calculation method can run with high speedup ration and parallel efficiency. The research will further advance the technology of anisotropic seismic exploration in coal fields, provide data for seismic inversion and give a theoretical support for coal mine disaster prediction.